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Abstract 


Seismic waves sensitive to the outermost part of the Earth’s liquid core seem to be 
affected by a stably stratified layer at the core-mantle boundary. Such a layer could have 
an observable signature in both long-term and short-term variations of the magnetic field 
of the Earth, which are used to probe the flow at the top of the core. Indeed, with the 
recent SWARM mission, it seems reasonable to be able to identify waves propagating in 
the core with period of several months, which may play an important role in the large-scale 
dynamics. 

In this paper, we characterize the influence of a stratified layer at the top of the core on 
deep quasi-geostrophic (Rossby) waves. We compute numerically the quasi-geostrophic 
eigenmodes of a rapidly rotating spherical shell, with a stably stratified layer near the 
outer boundary. Two simple models of stratification are taken into account, which are 
scaled with commonly adopted values of the Brunt-Vaisala frequency in the Earth’s core. 

In the absence of magnetic field, we find that both azimuthal wavelength and frequency 
of the eigenmodes control their penetration into the stratified layer: the higher the phase 
speed, the higher the permeability of the stratified layer to the wave motion. We also show 
that the theory developed by Takehiro and Lister 2001 for thermal convection extends 


to the whole family of Rossby waves in the core. Adding a magnetic field, the penetrative 
behaviour of the quasi-geostrophic modes (the so-called fast branch) is insensitive to the 
imposed magnetic field and only weakly sensitive to the precise shape of the stratification. 

Based on these results, the large-scale and high frequency modes (1 to 2 month periods) 
may be detectable in the geomagnetic data measured at the Earth’s surface, especially in 
the equatorial area where the modes can be trapped. 


1 Introduction 


The existence of a stably stratified layer at the top of the Earth’s liquid core can have conse¬ 
quences on different aspects of the Earth’s core dynamics. Such a layer can have a chemical 
origin with the accumulation of light elements coming from the growth of the inner core 
[Fearn and Loper 1981, Braginsky, 2007 Gubbins and Davies, 2013 or from chemical in¬ 


teractions with the mantle if the outer part of the core is under-saturated in oxygen and 
silicium |Buffett and Seagle, 20101. It can also have a thermal origin, if the heat flux at the 
core mantle-boundary is lower than the adiabatic heat flux in the core. With their recent 
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estimations of the thermal iron conductivity in the Earth’s core, de Koker et al. 2012 and 


Pozzo et al. [2012 claim that such a layer can be thermally stable at the top of the core. 

The strength of the stratification in such a layer is characterized by the Brunt-Vaisala fre¬ 
quency 


N = \ —--r~ 


9 dp 
p dr ’ 


(1) 


where p(r) is the density profile as a function of the radius, and g the acceleration of gravity, 
which is considered constant in the thin layer [e.g. Crossley, 1984 . A stably stratified layer 


corresponds to the case N 2 > 0, and a convective region to N 2 < 0. The comparison between 
N and the planetary rotation rate 11 allows to distinguish two stratification regimes: weak 
stratification when N < 212 and strong stratification when N > 212. 

The presence of such a stratified layer at the top of the outer core has been advocated 
by some seismologists to explain travel-time anomalies of seismic body waves sensitive to 
the outermost part of the Earth’s liquid core. It is however difficult to detect small density 
anomalies in this region, because body waves are also affected by the strong heterogeneities of 
the D” layer at the base of the mantle [Souriau , 2007|. Optim istic estimations give a thickness 
of the stratified layer between 50-100 km |Lay and Young 19901 and 300 km |Helffrich and 
Kaneshima, 2010 , with a weak stratification corresponding to 1V/212 ~ 0.1 in the latter case. 

; Coriolis and buoyancy force 

These inertia-gravity waves have been investigated 


In a stratified layer, waves in which the Coriolis and buoyancy forces play an important role 


1982 


for several decades in the geophysical and astrophysical contexts Olson 

1977 

Crossley and 

Rochester, 

1980, 

Crossley 

1984, 

Friedlander 

1985 

Dintrans et al. 1999 

. If we could detect 


them, they would also provide evidence for the existence of a stably stratified layer in the 
outermost part of the Earth’s core. The possible detection of these waves in the geodetic 
high frequency spectrum attracted attention some years ago, thanks to the development 
of superconducting gravimeters Melchior and Ducarme, 1986 Aldridge and Lumb 1987 


Unfortunately, the observations in the high frequency range do not require a stably stratified 
layer to explain the data Melchior et jffij 1988" , mainly because of the uncertainty on the 


frequencies of the waves [Rieutord 


1995 


The evolution of the magnetic field of the Earth 


also gives constraints on the stratification at the top of the core Gubbins, 2007 


Although still debated, and with characteristics such as depth and N not known, the 


et al. 

1997 

Lister and Buffett, 1998 Labrosse, 2003 . This layer may also modify the dynamics 

of the Earth’s 

liquid core. 

Zhang and Schubert 

1997L 

Takehiro and Lister 

2001 

2002 and 

more recently 

Nakagawa 

2011 studied the penetration of columnar convective motions of 


a rapidly rotating spherical shell in a stably stratified layer located at the outer boundary. 
Their results show that the stratification acts as an interface, which can be crossed by large- 
scale motions while small-scale ones are trapped below the layer when the stratification is 
strong enough. Such a layer has then the ability to partially hide geophysical flows inside the 
Earth’s core, which has important implications for core flow inversion. 

Indeed, using recordings of the magnetic field, it is possible to infer the flow at the core 
surface using inversion techniques Holme, 2007 . Assuming a quasi-geostrophic (columnar) 


flow, it is then possible to reconstruct the flow in the core interior [see e.g. Pais and Jault| 


2008J. Although there is evidence for quasi-geostrophic dynamics in the Earth’s core~[Gillet 


et al. 2011 , the presence of a stratified layer may prevent us to reconstruct part of the interior 
flow [Bloxham, 1990j. Moreover, a stratified layer can host its own eigemnodes, which may 
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Figure 1: (a) Geometry of the system, (b) Dimensionless Brunt-Vaisala frequency N(r)/£l 
as a function of the dimensionless radius r. The solid red and dashed blue curves correspond 
to model (13) and model (14) respectively, with Nq = fi. The thickness of the stratified layer 
is 140 km. 


then be distinct from the flow in the bulk of the core [Braginsky] 1998, 1999, Buffett, 2014 


In this study, we characterize the effect of a stratified layer on otherwise quasi-geostrophic 
inertial modes in the Earth’s core, also known as Rossby modes. Their time-scale is typically 
of a few months, and their signature may be actually seen in the geomagnetic data from the 
new SWARM mission Olsen et ah, 2013 . Using a numerical eigennrode solver, we compute 


the linear hydromagnetic modes in a simplified model of the Earth’s core including a stratified 
layer at the outer boundary. The effect of an imposed magnetic field on the modes is briefly 
considered in the case of a simple toroidal magnetic field depending only on the cylindrical 
radius s [Malkus. 1967 . 

The paper is divided as follows. Section [2] presents the mathematical modelling and 
section [3] the numerical methods used to solve the problem. The results are shown in section 
[4] with an emphasis on the shape of the eigenmodes and the dependence of their penetration 
length on governing parameters. Section [5] discusses our findings. We end the paper with a 
conclusion. 


2 Mathematical modelling 

2.1 Basic equations 

The geometry of the system is illustrated in the Figure [pi. We consider a rapidly rotating 
spherical shell of inner radius Ri and outer radius R 0 , surrounded by perfectly insulating 
mantle and inner core. We assume an aspect ratio 7 = Ri/R 0 = 0.35, similar to that of the 
Earth’s core. The shell, the mantle and the inner sphere are rotating at the same angular 
velocity Ct = 9z aligned with the unitary vertical axis z in cylindrical polar coordinates 
(s, (j>, z). We work in the rotating reference frame attached to the mantle, in cylindrical polar 
coordinates or in spherical polar coordinates ( r,0,q !>). The velocity v, the magnetic field B, 
the temperature T and the density p are subjected to small perturbations around a basic 
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state (Uq, Bq, Tq, po). Let us write 


( 2 ) 

(3) 

(4) 


[v, B] (r, t) = [U 0 , B 0 ] (r) + [u, b] (r, t), 
T{r,t) = T 0 (r) + 0(r,t), 
p(r,t) = po (1 - a©(r,t)), 


where a is the volume coefficient of thermal expansion and (u, b, 0 , —poaO) the perturbations 
of small amplitude. For the sake of simplicity, we impose Uo(r) = 0, which corresponds to 
a solid-body rotation in the inertial frame. We also assume that the background density p o 
is constant, neglecting the variations due to adiabatic compression across the shell Rieutord 
19951. We take into account the density perturbation only in the gravity term with the 


the 


Boussinesq approximation, considering a gravity field of the form g = — g Q r r with g 0 R 0 
acceleration of gravity at the outer boundary and r the unitary radial vector. 

To keep the model simple, we suppose that the background magnetic field Bo in the shell 
is either an axial field 

B 0 = Bq z, (5) 


or the Malkus field [Malkus 1967] 


The latter is a purely toroidal and azimuthal held which increases with the cylindrical radius 
s = rsin(9. This assumption is consistent with our choice for the background velocity, since 
no basic velocity Uo is needed for the basic state to be in equilibrium (the resulting Lorentz 
force is balanced by the pressure gradient). 

We use R 0 as the length scale, fl _1 as the time scale, pqR 0 VL 2 as the pressure scale, 
Ll 2 R 0 /ag 0 as the temperature scale and y/poPo QRo as the magnetic scale. As we are interested 
by the eigenmodes, we neglect the non-linear terms and the dimensionless equations satisfied 
by the perturbations are 


<9u 

dt 


+ 2 z x u = —VP + Le [(V x b) x bo + (V x bo) xb] + r0 r + E V 2 u 


<90 N 2 (i 

~dt = 

dh 

— = Le V x (u x bo) + E v V 2 b. 


u-f + J^V 2 0, 
Pr 


(7) 

( 8 ) 

(9) 


Five governing parameters appear in the previous equations, namely the Ekman number E, 
the magnetic Ekman number P, ; , the Prandtl number Pr, the Lehnert number Le and the 
dimensional Brunt-Vaisala frequency N. They are defined by 


E = 


nR 2 ' 


T? _ '< 

71 ~ QR 2 ’ 


Pr = -, 

K 


Le = 


Bn 


\JPoPo LlR 0 


and N(r) = \/ag 0 (10) 


with v, k and g respectively the viscous, thermal and magnetic diffusivities and Bq the 
strength of the background magnetic field. In the computations, we take Pr = 0.1 — 10 and 
Le = 10 -4 . Because of numerical constraints, we choose E = E v = 10 - '. 
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Equations Q, ([8]) and ([£]) are completed with boundary conditions on the velocity, the 
temperature and the magnetic field. We impose on both shells the no-slip boundary condition 
for the velocity field 

u = 0 (11) 


and a constant heat flux 


30 „ 

Ih ~ ° * 


( 12 ) 


The magnetic field boundary conditions are the continuity of b and the continuity of the 
radial component of V x b. The spherical harmonic expansion (introduced below) allows to 
write this as a boundary condition. 


2.2 Basic state for the stratification 


The Brunt-Vaisala frequency N(r ) in the core is highly speculative. Following Lister and 
Buffett 11998 , let us assume that the spherical shell is divided in two parts (Figure [TJi). A 
stably stratified layer with N 2 (r) > 0 is located above a convective core of dimensionless 
radius d, which is supposed to be neutral with respect to the stratification (N 2 (r) = 0). This 
assumption filters out the coupling between waves and convection |Rieutord, 19951. We set 
the radius d = 0.96, which corresponds to a stratified layer of thickness 140 -150 km. This 
thickness belongs to the range estimated by seismic studies 


Hclffrich and Kaneshima 2010 


It is also inferred by Buffett 2014 to match the geomagnetic secular variation of the axial 


dipole with slow Magneto-Archimedes-Coriolis (MAC) waves. Following Takehiro and Lister 
2001, 2002 and Nakagawa 2011 , we adopt 


rr 

{ r — d\ 

i ~ ~ 

1 — tanh ( - ) 

{ 2 

L V a J\ 


(13) 


which is illustrated in Figure [TJx Nq is the Brunt-Vaisala frequency at the outer boundary 
and a = 0.005 the transition thickness between the neutral core and the stratified layer. N(r) 
is thus a smooth profile in the core with a roughly constant value N(r) ~ No in the outer 
layer. Different values of No are taken into account. In addition, we also consider a profile 
similar to the one used by Rieutord 119951 and Buffett] [[2014; 


i.e. 


N(r ) = 


2No{r — d)/(l — d) V r > d 
0 V r < d 


(14) 


It is linear in the stratified layer, with a mean value Nq and a maximum value 2Nq at the 
outer boundary. 


3 Numerical modelling 

3.1 Spectral method 

As they are solenoidal, u and b are decomposed into poloidal and toroidal parts 

u = V x V X (Ur) + V x (Vr), 
b = V x V x (Gr) + Vx (Hr), 
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(15) 

(16) 



































with (U, G ) and (V, H) the poloidal and toroidal scalars. The scalars (U, V, G, H , 0) are then 
expanded on orthonormalized spherical harmonics of degree l and order rn, as 

OO 

X] [ur, vr , ©n w ytw, 0) exp^, a?) 

l=m 
00 

£ [GT- «T] M V(«, 0) exp(At). (18) 

l=m 

Since the basic state (Uo,B 0 ) is axisymmetric, there is no coupling between Fourier modes 
m in the equations 0 to 0- Thus, one can seek solutions (u, 0, b) with a single azimuthal 
wavenumber m. As we are interested in the eigenmodes of the system, we assume that the 
time dependence of all perturbations is proportional to exp(At), with the complex eigenvalue 


[U, V, 0] (r ,t) = 
[G,H,](r,t) = 


A = t + iu) (19) 

where r is the growth rate of the mode and u its frequency, both non-dimensional. 

Thanks to the spectral decomposition, the boundary conditions are easy to express. On 
the two boundaries, the no-slip condition becomes 


U, 


m 


d U™ 
dr 


= vr = 0 


and the constant heat flux condition 


d ©m 

dr 


= 0 . 


( 20 ) 


( 21 ) 


With the insulating boundary condition, the magnetic field matches a potential field in both 
the mantle and the inner sphere, resulting to 


7 

—1 - GV 1 = HJ n = 0 (22) 

dr r 

on the inner boundary (r = 7 ) and 

Arim I 1 1 

+ ——G™ = HJ 71 = 0 (23) 

dr r 1 ‘ 


on the outer boundary (r = 1 ). 


3.2 Symmetry considerations 

The equations 0-0 are symmetric or anti-symmetric with respect to the equatorial plane. 
This implies that solutions can be characterized by their symmetry with respect to the equator 
[Gubbins and Zhang, 1993 . 

In the following, we focus on equatorially symmetric u and 0, because quasi-geostrophic 
flows have this equatorial symmetry. They satisfy 

[u r ,ug,Uti,](r,9,(f>) = [u r , -itg, u^](r, 7T - 9,(f>), 

©(r,M) = 0 (r, 7 r — 6, cf>). 
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(24) 

(25) 












Figure 2: Sparse matrices A and B in the non-magnetic case. L max = 6 and N r = 20. The 
black lines separate the poloidal and toroidal components of the velocity and the temperature 
field. The first line of blocks represents the evolution equation of the poloidal component of 
the velocity (coupled to the toroidal flow by the Coriolis force and to the temperature by the 
Buoyancy force), the second line is the equation for the toroidal velocity, and the last line is 
the equation for the temperature (coupled to the poloidal flow only in the stratified layer). 


Taking the symmetry of Yj m into account, the spherical harmonic expansions for the poloidal 
and toroidal components of u (and b) alternate to match the equatorial symmetry. Note also 
that the Coriolis force couples U™ and V)"']. Following 
u and 0 is given by 

(TT m V m £■)’ 

\ u l ) v l+ii 


Schmitt 


2010 


the symmetric part of 


(26) 


for a given azimuthal number m and with l = m,m + 2,m + 4,... 

Because u is symmetric, the symmetry of b is the same as the one of Bo: b is equatorially 
symmetric for the Malkus field, while it is equatorially antisymmetric for the axial magnetic 
field. Hence, for the magnetic field we add to the state vector 


or ( Gr +1 ,HD 


respectively for symmetric or antisymmetric imposed fields. 


(27) 


3.3 Numerical implementation 


Substituting from (17) and (18) into the governing equations ([7]) to leads to a set of 
linear differential equations which forms a generalized eigenvalue problem with its boundary 
conditions. To get numerical solutions, the spherical harmonic expansions are truncated at 
the harmonic degree Lm and the differential operators are represented by a second order 
finite difference scheme on an irregular mesh of N r levels, adjusted to properly resolve the 
Ekrnan boundary layers. The eigenvalue problem is now in a matrix form as a complex, 
non-hermitian and generalized eigenvalue problem 


HX = A BX. 

A is the eigenvalue and X is the associated symmetric eigenvector given by 

X=(7/r(r i ),^i(^),©r^),Gr(r i ),^ 1 (r0) t , 


(28) 


(29) 
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with i = 0,1,..., N r and l = m, m + 2, m + 4,... for a symmetric eigenmode permeated by 
a symmetric imposed field. The two matrices A and B , collecting the terms on the left-hand 
side and right-hand side of the equations, have approximately (5 Lm-^t/2) 2 elements. In order 
to have good spatial and spectral resolutions, we use Lm ~ 350 and N r — 850 for Ekman 
number E = 1CP 7 . Matrices A and B have a large number of elements, but most of them are 
zero. They are stored as sparse matrices in order to reduce the required memory (see Figure 


§• 


Following Rieutord et al. 1997 and Dintrans et al. [1999 , the linear eigenvalue prob 


lem (28) is solved using the shift-and-invert spectral transformation: instead of solving the 
problem (28), the related problem 


(A - crR) _1 RX = /iX 


(30) 


is solved with a the shift and the new eigenvalue /i related to A by 

1 

l 1 = 


A — a 


(31) 


With this spectral transformation, it is easier to get the least damped quasi-geostrophic modes, 
imposing a = iuj with 1. Only eigenvalues with absolute residuals |(A — aB)~ 1 B~K — 
MX| < 10 12 are accepted as good eigenvalues. 

We wrote a Python code, nicknamed Singe (Spherical INertia-Gravity Eigenmodes), that 
relies on the SLEPc library [Hernan dez et al., 2005] to solve the eigenvalue problem using the 


restarted Krylov-Schur method |Hernandez et al., 2009 together with the shift-and-invert 
strategy. SLEPc is an open-source scientific library developed to efficiently solve large and 


sparse eigenvalue problems on parallel computers. It is itself based on PETSc |Balay et al. 


1997 2014 


another open-source library intended to solve scientific problems modelled by 
Finally, Singe uses the SHTns library 


partial differential equations 
spherical harmonic transforms. 


Schaeffer, 2013 for the 


The Singe code was benchmarked with the theoretical solutions of Zhang et al. 12001 


Liao et al. 2001 and Zhang et al 


shell geometry, we used the numerical results of Rieutord et al. 1997 


2003 in the full sphere. As benchmarks for the spherical 

without stratification 


and of Dintrans et al. 1999 with a constant stratification in the shell. 


Singe is available at https://bitbucket.org/nschaeff/singe as free software. 


4 Results 


4.1 Description of eigenmodes 

To study the behaviour of quasi-geostrophic inertial modes in the presence of a stratified layer, 
we first focus on the non-magnetic case. Characteristics of some modes are gathered in Table 
|T] Each mode is specified by a triplet ( u>,m,n c ), with m the azimuthal wave number and 
n c the number of zero-crossings of a velocity component ( u s or u^) along a cylindrical radial 
profile (measuring latitudinal variations). The dimensionless frequency cj, with correspond¬ 
ing periods ranging from one to several months, is always negative because the eigenmodes 
propagate eastwards. It is a characteristic of Rossby waves in the Earth’s core geometry [e.g. 
Busse, 1970|. Indeed, looking at the different forces, we find that the modes are dominated 
by the Coriolis force, with some exceptions in the stably stratified layer (see below). The 
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LV 


Tl c 

O 

II 

£ 

Cl 

II 

S? 

II 

to 

No = io n 

solid SF 

solid NS 

8 

-0.0046 

-0.0051 

-0.0054 

-0.0060 

-0.0051 

-0.0043 

3 

-0.014 

-0.017 

-0.020 

-0.023 

-0.015 

-0.014 

13 

-0.0039 

-0.0044 

-0.0048 

-0.0051 

-0.0042 

-0.0036 

7 

-0.019 

-0.020 

-0.023 

-0.025 

-0.020 

-0.019 

13 

-0.0069 

-0.0075 

-0.0080 

-0.0084 

-0.0072 

-0.0066 

8 

-0.021 

-0.024 

-0.029 

-0.031 

-0.022 

-0.021 

14 

-0.0079 

-0.0084 

-0.0093 

-0.0093 

-0.0083 

-0.0076 

10 

-0.018 

-0.020 

-0.022 

-0.024 

-0.018 

-0.017 


Table 1: Wave number m (azimutal), approximative number of quasi-geostrophic columns 
along a meridian n c , dimensionless frequency u and period in months of some numerical 
eigenmodes for four values of No, without magnetic field. Computations at E = 10~ 7 and 
Pr = 0.1. IV(r) is given by (13). In addition, we have also computed the same modes by 
replacing the stratified layer by a solid layer, with stress-free (solid SF) or no-slip (solid NS) 
boundary conditions. 


Coriolis force tends to align the flow along the spin axis to create Taylor columns, and the 
modes are mainly quasi-geostrophic inertial modes (almost invariant along the rotation axis). 
The frequency depends also on the length scale of the modes. For a given azimuthal number 
m, high frequency modes are large scale modes with a small number of columns n c , while low 
frequency are small scale modes with a large number of columns. 

The precise shape of the modes depends on the stratification, which may influence the 


1997, Takehiro and Lister 


penetration of the flow in the outer layer Zhang and Schubert 
2001 . However, the overall pattern of the modes in the neutral part of the core smoothly 


evolves with the stratification strength Nq. For commonly adopted values of the Brunt- 
Vaisala frequency Nq, the number of columns of a mode remains approximatively the same, 
while the exact position and thickness of the columns change near the outer layer. A smooth 
increase of IVo allows us to study the evolution of the modes. It can be done with an iterative 
process to follow a given mode defined by (m, n c ). As observed in Table [lj uj slowly increases 
with No- This slow increase of the frequency uo cannot be attributed to a geometrical effect. 
Indeed, when replacing the stratified layer by a solid layer, the variation in frequency is much 
smaller, especially at large wave-numbers (see Table [I]). 

The cylindrical components u s and u^ in the equatorial plane of some modes are shown 
in Figure [3j It is worth noting that lower frequency modes are spiralling more than higher 
frequency modes. The spiralling appears because the viscosity is overestimated in our com¬ 
putations (see appendix [A| for more details). Finally, the kinetic energy is dominated by the 
azimuthal component of the velocity u which is more concentrated in the equatorial region 
when the frequency is higher. This phenomenon, called trapping, is known for inertial waves 
Zhang, 1993 as well as for inertia-gravity waves [Crossley 1984 Friedlander, 1985 


E < 10“^) we observe that this trapping becomes more efficient. Unfortunately, it is then 
expensive to resolve numerically, because it needs enough grid points to capture the boundary 
layer and enough spherical harmonics to describe the trapping. For these reasons, we set the 
Ekrnan number to 1CU 7 in our computations. 

Let us now switch on the background magnetic field Bo given by expression (J6|) . Without 


When 
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Figure 3: Contours of the cylindrical components u s and u<f, of four numerical eigenmodes 
from Table [I] (without magnetic field) in the plane (m, |w|). lo is the dimensionless frequency. 
Contours are shown in the equatorial plane as seen from above. The vertical component u z 
is not shown because it vanishes in the equatorial plane. The velocity is normalized by the 
maximum value of u s . The black dashed circles represent the transition between the neutral 
core and the stably stratified layer. The modes were computed at E = 10 ' and Pr = 0.1. 
N(r) is given by (13) and Nq = 12. (a) High-frequency mode m = 1. (b) Low-frequency mode 
m = 1. (c) High-frequency mode m = 6. (d) Low-frequency mode m = 6. 
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stratification, a background magnetic field splits the modes into ’’slow” and ’’fast” magneto- 
Coriolis (MC) waves, the latter being essentially pure inertial waves [see e.g. 

2003]. 

We computed these fast modes with Pr = 0.1 and a Lehnert number Le = 10 -4 corre 


Zhang et al. 


sponding to Bo = 3mT [in the range estimated for the Earth’s core, see Gillet et al. 
For the Malkus field (eq. [6” 


2010 


we were able to compute the modes down to E = E v = 10“ 
Computation of the modes with an axial field (eq. [5]) was numerically more demanding, 
probably because of the additional coupling terms in the corresponding matrix. Indeed, in 
that case we were able to compute the modes down to E = E v = 10 -5 only. For azimuthal 
wave-numbers ranging from m = 1 to m = 26, the frequency cj changes by less than 1% 
and the shape of the modes is indistinguishable from the ones obtained without magnetic 
field. Note however that the overall damping rate of the mode is affected, especially when 


10 5 , but we expect it to remain small in the Earth’s core, where E^ ~ 10 


-10 


E = E v = 

This holds for every eigenmode considered here, for all values of the Brunt-Vaisala frequency 
No investigated here, and for the two magnetic field models (Malkus and axial). It shows 
that the eigenmodes considered in this study are barely modified by the background magnetic 
field, because they are mainly ’’fast” quasi-geostrophic inertial modes. These results are in 
agreement with the theoretical predictions of Malkus 1967 and Zhang et al. [2003 for a 
homogeneous fluid contained in a full sphere. Furthermore, Canet et al. [2014 considered a 
quasi-geostrophic model of magneto-Coriolis (MC) waves and showed that various shapes of 
a toroidal background magnetic fields do not affect the fast modes either. 

Since our quasi-geostrophic modes are not affected significantly by the magnetic field, with 
or without a stratified layer, we will no longer consider the magnetic field in the remainder 
of this paper, effectively solving only equations 0 and ([8]) with Le = 0. This simplifies the 
numerical study quite a bit. Note that although the fast branch considered here is unaffected, 
the slow branch will most certainly be affected by the magnetic field, but this is beyond the 
scope of this paper. 


4.2 Penetration into the stratified layer 

In order to quantify the penetration of the quasi-geostrophic modes into the stratified layer, 
we compute profiles at a given distance s from the rotation axis. Figure [4] shows the profiles 
of the azimuthal average of the quasi-geostrophic kinetic energy u 2 + u ^ along the axial 
direction, for three values of No and two wave numbers m. For a given Nq, the modes have 
approximately the same frequency. The points are chosen along the column located at the 
cylindrical radius s = 0.7, where the energy is non-zero in the equatorial plane. The sharp 
bump before the decrease of the kinetic energy, observed on all curves, corresponds to the 
Ekman layer. When No = 0.1 P, the Taylor column penetrates into the outer layer. The flow 
just below the Ekman layer is nearly the same as that in the neutral core. When No = P, the 
kinetic energy decreases in the layer, although it does not vanish before the Ekman layer is 
reached. When No = 2 P, the kinetic energy is further reduced, and even more so for higher 
azimuthal modes m. 

From these profiles (Fig. [4]), we can define a transmission coefficient T. It is our key 
parameter to compare quantitatively the penetration of the flow inside the stably stratified 
outer layer, for eigenmodes of various shapes and frequencies. Let us define T as the ratio of 
the value of the square root of the azimuthal average of the quasi-geostrophic kinetic energy 
u 2 s + u^, taken just under the Ekman boundary layer (Figure®, to that in the equatorial 
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m = 1 m = 8 



Figure 4: Distribution of the azimuthal average of the quasi-geostrophic kinetic energy u^+ui, 
along a quasi-geostrophic column located near the cylindrical radius s = 0.7, for some high 
frequency modes of the Table [lj computed at E = 1CP 7 and Pr = 0.1. N(r) is given by eq. 
(13). The horizontal axis z — zi ayer represents the difference between the vertical coordinte z 
along the column and the coordinate of the bottom of the stratified layer zi ayer . The vertical 
dashed line marks the position of the transition between the neutral core and the stably 
stratified layer. For each m, the modes have approximately the same shape in the neutral 
core. The bumps in the energy at the right-hand side of the curves show the effect of the 
Ekrnan pumping. The curves are normalized by the value of their energy in the equatorial 
plane (not shown). 
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Figure 5: Maps of the amplitude of the transmission coefficient T as a function of the az¬ 
imuthal number m and the dimensionless frequency |w| of the modes. Top row: layer with 
constant Brunt-Vaisala frequency (eq. 13); bottom row: layer with linear Brunt-Vaisala fre¬ 
quency (eq. 14). The Brunt-Vaisala frequency is set to Nq = If on the left column and 
Nq = 2on the right column. Values are computed at the cylindrical radius s = 0.7. T 
is interpolated between the black dots which represent numerical modes for which there are 
quasi-geostrophic columns at s where the kinetic energy is non-zero. In the two models, the 
stratified layer has the same thickness (140 - 150 km). 
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Figure 6: Dependence of transmission coefficient T with the dimensionless phase velocity 
= \u\/m of the modes, for various cylindrical radius s and Brunt-Vaisala frequency Nq. 
Computations were performed at E = ICC 7 , Pr = 0.1 and Nq = D. T is computed automat¬ 
ically for all the modes. Only the results for constant Brunt-Vaisala frequency (eq. 
shown. 


13) are 


plane. Figure [5] shows the distribution of transmission coefficients T defined along the quasi- 
geostrophic column at a given cylindrical radius (s = 0.7). Modes of different m and different 
frequencies u are represented. For a given Nq, large-scale modes are less affected than small- 
scale modes which cannot penetrate deep into the layer. This behaviour is reinforced when 
Nq increases, showing that the stratification clearly acts as a low-pass filter: large-scale 
eigenmodes (high frequencies, low m) penetrate deeper inside the outer layer than small-scale 
eigenmodes (low frequencies, large m). 

Because N(r ) is not known in the Earth’s core, Figure [5] also compares the effect of two 
stratification profiles (13) and (14). The modes can penetrate slightly deeper in the layer 
with the linear profile, but the overall distribution is similar. Therefore, the exact shape of 
N(r ) does not seem to affect the leading order behaviour for a layer with a given thickness. 
Finally, we have also observed that the thicker the layer, the less the modes penetrate into 
the outer layer, which is the expected behaviour |Takehiro and Lister, 20011. 

The lines of constant transmission coefficient T displayed in Figure [5] suggest that a single 
parameter may control T. Figure [6] shows the transmission coefficient T as a function of the 
dimensionless angular phase velocity = \oj\/m for profiles taken on modes with various m 
and cj, three cylindrical radii s, and two values of Nq. The modes with different m and uj 
collapse quite well, although a residual dependence on s (the location of the quasi-geostrophic 
column in the shell) and Nq is now apparent. Columns at high and mid latitudes can penetrate 
slightly better in the outer layer than columns at low latitudes. At a given s , modes with a 
larger phase velocity are in general less affected by the outer layer than the others. There 
are also some departures from the general trend, which are likely to come from the automatic 
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n 0 = n N 0 = 2n 

u = -0.0115 w = -0.0138 

Figure 7: Distribution of the azimuthal component u^ of the velocity held for two large-scale 
eigenmodes m = 1 in a meridional plane. The dashed line represents the bottom of the 
stratified layer. The amplitude is arbitrary. Computations at E = 10”' and Pr = 0.1. 


procedure to follow modes as the parameters are changed and the automatic evaluation of T. 


4.3 Equatorial trapping 


A peculiar behaviour is observed near the equator. There, we observe that modes of di¬ 
mensionless frequency |w| > 0.02 and of small wave numbers m < 3 can penetrate inside 
the equatorial region, where the azimuthal component u$ can be equatorially trapped in the 
stably stratified layer. This phenomenon is illustrated in Figure [7j In this case, the modes 
still have uo < 0, reminding of Rossby waves, but they are now governed by both Coriolis 
and buoyancy forces. This trapping has been predicted by theoretical studies Friedlander 


and Siegmann, 1982, Crossley, 1984|. Note however that a strong dipolar magnetic field may 


release the trapped modes Bergman, 1993 


4.4 Thermal vs chemical stratification 

Thermal and chemical buoyancy are governed by the same advection-diffusion equations, 
the main difference being the very different value of their respective diffusivity. Estimates 
of diffusivities in the Earth’s core suggest that representative values of Prandtl number are 
Pr >0.1 for thermal diffusivity and Pr > 10 for the diffusivity of light elements in the 
core. We can thus study the effect of thermal or chemical stratification by adjusting the 
value of the Prandtl number in our model. Figure [8] shows the dependence of the absolute 
dimensionless frequency |cu| and of the transmission coefficient T as a function of Pr, for some 
high frequency modes of different wave numbers m. For Pr > 10, an asymptotic regime is 
reached where both ui and T do not depend on Pr any more. In any case, Pr has only a 
small effect on the penetration of the columns in the outer layer. 
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Figure 8: Evolutions of the absolute dimensionless frequency |w| and of the transmission 
coefficient T as a function of the Prandtl number Pr of some quasi-geostrophic modes. Com¬ 
putations at fixed E = 10 - ' and for model (13) with Nq = fi. 


5 Discussion 

5.1 Comparison with previous studies 

We have shown that a sufficiently stratified layer can affect the penetration of quasi-geostrophic 
inertial modes based on their phase speed. A theoretical value for the penetration distance 5 
has been derived by Takehiro and Lister 2001 , using a local plane layer approximation, 


5th = 


1 


20 


N 0 y/k 2 + l 2 


(32) 


where k and l are the wavenumbers in the azimuthal and radial directions respectively. 
Because l is difficult to estimate, they further simplify their expression assuming isotropy 
l = k = m/s for a global neutral convection mode of azimuthal wavenumber m. Although 
this hypothesis is questionable (and we will discuss a better one below), they find good agree¬ 
ment between the penetration length obtained in numerical simulations and their estimate 

_ 20 
0 m. — T7“ 


No \J 2 to ’ (33) 

as long as b < 0.4. However, this last expression 5 m translates into a transmission coefficient 
T that depends only on the azimuthal wavenumber m, in contradiction with our findings that 
also exhibit a strong dependence with the frequency u (see Fig. [5]). 

We have found it difficult to measure reliably a penetration distance 5 on the profiles 
shown on Figure [4j where the amplitude seems to decay exponentially in the lower region of 
the layer and then saturate to a constant value near the solid boundary. We therefore defined 
a transmission coefficient, which has the drawback of being dependent of the layer thickness, 
but the advantage of being well defined, even for arbitrarily shaped profiles. 

The modes we have studied are very close to global Rossby modes, which have the following 


dispersion relation in a local plane approximation without stratified layer Greenspan 1968 : 
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(34) 


cj = 2 np 


k 


k 2 +1 2 


where /3 = — s/(l — s 2 ) for a full sphere, and all length scales are normalized by the radius of 
the Earth’s core Rq. Ignoring the small variations of u with N (see table [I]), we can use this 
expression to estimate k 2 + 1 2 from u and k = m/s. Equation (32) then becomes 


SRossby ~ ^ \j 212(1 S 2 ) 'J - ^ } 


■\J 2Q Ic^l, 


(35) 


where 6 is the colatitude and c^ the phase speed in the azimuthal direction. As shown on 
Figure [9] this last expression for 6 that assumes the local dispersion relation of Rossby waves 
(eq. [34] ) seems to control all the variability in the transmission coefficient T. 

We have thus shown that the local plane-layer theory of Takehiro and Lister 12001 for 
thermal convection extends to the whole family of Rossby waves (quasi-geostrophic inertial 
waves) in the core. This implies a transmission coefficient T for Rossby waves that depend 
only on phase speed = oj/rn and on colatitude 6 in addition to the control parameters Nq 
and Q. 


5.2 Possible observations in the geomagnetic data 

It is worth noting that the the absolute frequency of the waves increases with Nq (Table 
[I]). This sensitivity of eigenmode frequencies to the Brunt-Vaisala frequency may be a way 
to probe for the existence of a stably stratified layer below the core-mantle boundary, but 
other details may also alter the frequency of the modes (ellipticity, mean flow, ...), and it will 
probably be difficult to tell them apart. 

In addition, the relatively high frequencies of our modes are very difficult to extract from 
the residual noise in the current geomagnetic models. However, the ability of the three 
SWARM satellites to better separate external and internal fields Olsen et ah, 2013 , may 
enable us to detect the signature of these modes in the geomagnetic data. 

We have computed the secular variation induced by our modes acting on a model of the 
present geomagnetic field Finlay et ah, 20101. The magnetic signatures of large-scale (m < 3) 
and high frequency modes (periods shorter than two months) concentrates in the equatorial 
area, where they are amplified (see Figure [8]). If these waves are detected, the width of the 
equatorial belt in which they appear will also unveil the depth of the stratified layer. 


5.3 Influence of the magnetic field 

The quasi-geostrophic fast modes are only very weakly sensitive to the simple magnetic fields 
considered here (the toroidal Malkus field and the uniform axial field), for magnetic field 
strength in the estimated range for the Earth’s core {Bq = 3mT). In particular, the magnetic 
field does not change their ability to penetrate the stratified layer. It comes from the fact that 
these modes evolve on a time scale much faster than that of the Alfven waves, as measured by 
the small Lehnert number of the Earth (Le = ICE 4 ). Note that the Lehnert number depends 
on a length-scale t. Le{t) = Bq/ y/poHo fi£, so that at small enough length-scales the magnetic 
field will affect even the fast modes. For the Earth however, this length-scale is beyond what 
can be observed | Gillet et al. 
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Figure 9: Transmission coefficient T as a function of theoretical penetration depth 5r oss b y 
(equation |35[), normalized by the stratified layer thickness 1 — d. Top: constant stratification 


profile (eq. [l3| ). Bottom: linear stratification profile (eq. 14). All transmission coefficients 
obtained for various position (s = 0.5, s = 0.65, s = 0.75) and two stratification strength 
(N = and N = 2 Cl) collapse quite well. 
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We expect that the fast Rossby modes remain insensitive to more complex magnetic 
fields. Unfortunately, our numerical code Singe is not well suited to take into account complex 
imposed fields, because they lead to an increasing number of coupled harmonic degrees, which 
are cumbersome to express and translate into matrix coefficients. A possible solution would 
be to use computer algebra systems (CAS) to write down the coupling operators between 
spherical harmonics in a symbolic way [as in Schmitt 2010), and then to export them into 


a language suitable for numerical computations. Another more flexible approach is to use 
so-called matrix-free methods, which avoid to explicitly form the matrices, but require to 
provide appropriate preconditioners. These methods may come forward in the forthcoming 
years. 


6 Conclusion 


We have studied the possible penetration of quasi-geostrophic modes in a stratified layer at the 
top of the Earth’s core. The ability to penetrate depends on both the frequency |cu| and the 
azimuthal wave-number m: faster modes penetrate more easily than slower modes. Combining 


the finding of Takehiro and Lister 2001 and the dispersion relation of Rossby modes, we put 


forward a single parameter that accounts for all the variability in the penetration capability 
of the modes (eq. [35] ) . 

For an Earth-like magnetic field strength, the two simple field geometries we have consid¬ 
ered here do not affect the shape or the frequency of our fast Rossby modes in any significant 
way. This is due to the small Lehnert number Le = 10~ 4 which ensures that inertial waves are 
much faster than Alfven waves at large spatial scales. Although more complex magnetic field 
should be considered, it is tempting to conjecture that these fast Rossby modes are immune 
to any magnetic held. 

Large scale modes at periods of a few months are slow enough for their magnetic signature 
to be detectable in principle, and can penetrate a stratified layer with Nq = 2f2. Close to the 
surface of the core, the most intense motions associated with the modes in the presence of a 
stratified layer are located close to the equator, where they seem to couple with a trapped 
mode (see Figure[8]). This is where a magnetic signature may be generated and detected. Such 
a localized signature would provide observational evidence and constraints on the hypothetic 
stratified layer. 


Acknowledgments The authors would like to thank the whole geodynamo team for fruit¬ 
ful discussions. The most demanding computations were performed using the Froggy plat¬ 
form of the CIMENT infrastructure (https : //ciment .ujf-grenoble . fr), supported by the 
Rhone-Alpes region (GRANT CPER07_13 CIRA), the OSUG@2020 labex (reference ANR10 
LABX56) and the Equip@Meso project (reference ANR-10-EQPX-29-01). Plots were pro¬ 
duced using matplotlib [Droettboom et al., 2014, http://matplotlib.org]. ISTerre is part 


of Labex OSUG@2020 (ANR10 LABX56). 


A Effect of viscosity on quasi-geostrophic inertial modes 

Our study includes a small viscosity, which results in a non-zero damping of the modes, a 
small effect on their frequency, but also a significant effect on their shape, which is spiralling 
in the presence of sufficient viscosity. 
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Figure 10: Frequency and damping rate of the quasi-geostrophic pure inertial mode 
m = 3, as a function of the Ekrnan number. The green dashed line correspond to an 
scaling, the expected asymptotic damping due to the Ekrnan layers in a full sphere 
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As shown in Figure [TO], the asymptotic frequency and damping rate are approached only 

shows that the spiralling shape vanishes at the 


11 


when E ~ 10 -9 . Coincidentally, Figure 
same very low Ekrnan number. 

This contrasts with thermal convection in a rapidly rotating sphere, where tilted or spi¬ 
ralling structures are a feature of the onset [e.g. Dormy et al. 2004 Takehiro, 2008, and 


references therein]. The main reason being that, at the scale at which the convection first 
occurs, the viscosity is never negligible. 

Note also that this does not depend on the boundary conditions (stress-free or no-slip), 
showing that the spirallization is an effect of the viscosity in the bulk. 

All this shows that very low Ekrnan numbers are needed to have an asymptotic structure 
of quasi-geostrophic inertial waves. 
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